/****************************************************************************
 * Copyright (C) 2009-2010 SciTouch LLC
 * 
 * This file is part of Indigo toolkit.
 * 
 * This file may be distributed and/or modified under the terms of the
 * GNU General Public License version 3 as published by the Free Software
 * Foundation and appearing in the file LICENSE.GPL included in the
 * packaging of this file.
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 ***************************************************************************/

#include "base_cpp/auto_ptr.h"
#include "math/algebra.h"
#include "base_cpp/array.h"
#include "graph/graph_highlighting.h"
#include "graph/graph_decomposer.h"
#include "graph/filter.h"
#include "graph/graph_affine_matcher.h"
#include "graph/edge_rotation_matcher.h"
#include "molecule/molecule_substructure_matcher.h"
#include "molecule/molecule.h"
#include "molecule/molecule_arom_match.h"
#include "molecule/molecule_stereocenters.h"
#include "molecule/molecule_3d_constraints.h"
#include "molecule/molecule_neighbourhood_counters.h"
#include "molecule/molecule_pi_systems_matcher.h"
#include "graph/graph.h"

MoleculeSubstructureMatcher::MarkushContext::MarkushContext (const Molecule &query_, const Molecule &target_) : 
TL_CP_GET(query),
TL_CP_GET(query_marking),
TL_CP_GET(sites),
depth(0)
{
   int i;

   query.cloneMolecule(query_, 0, 0);

   const MoleculeRGroups &rgroups = *query.getRGroups();
   sites.clear();

   for (i = rgroups.begin(); i < rgroups.end(); i = rgroups.next(i))
      sites.push(rgroups.atomIdx(i));

   query_marking.clear_resize(query.vertexEnd());

   for (i = query.vertexBegin(); i < query.vertexEnd(); i = query.vertexNext(i))
      query_marking[i] = -1;
}


MoleculeSubstructureMatcher::MoleculeSubstructureMatcher (Molecule &target) :
_target(target),
TL_CP_GET(_3d_constrained_atoms),
TL_CP_GET(_unfolded_target_h),
TL_CP_GET(_query_h_counters),
TL_CP_GET(_target_h_counters)
{
   vertex_equivalence_handler = NULL;
   use_aromaticity_matcher = true;
   use_pi_systems_matcher = false;
   highlighting = 0;
   _query = 0;
   _am = 0;
   _pi_systems_matcher = 0;
   match_3d = 0;
   rms_threshold = 0;
   find_all_embeddings = false;
   
   _query_nei_counters = 0;
   _target_nei_counters = 0;

   calcTargetHCounters(_target, _target_h_counters);
   
   // won't ignore target hydrogens because query can contain
   // 3d features, hydrogen isotopes, etc.
}

bool MoleculeSubstructureMatcher::_shouldUnfoldTargetHydrogens (const Molecule &query)
{
   int i, j;

   for (i = query.vertexBegin(); i != query.vertexEnd(); i = query.vertexNext(i))
   {
      // lone hydrogen?
      if (query.getAtom(i).pureHydrogen() && query.getVertex(i).degree() == 0)
         return true;

      const QueryAtom &qa = query.getQueryAtom(i);

      // atom list that can include hydrogen?
      if (qa.type == QUERY_ATOM_LIST && qa.listContains(ELEM_H))
         return true;

      // atom not-list that can include hydrogen?
      if (qa.type == QUERY_ATOM_NOTLIST && !qa.listContains(ELEM_H))
         return true;
   }

   MoleculeRGroups *rgroups = query.getRGroups();

   if (rgroups != 0)
   {
      int n_rgroups = rgroups->getRGroupCount();
      for (i = 0; i < n_rgroups; i++)
         for (j = 0; j < rgroups->getRGroup(i).fragments.size(); j++)
            if (_shouldUnfoldTargetHydrogens(*rgroups->getRGroup(i).fragments[j]))
               return true;
   }

   return false;
}

void MoleculeSubstructureMatcher::calcQueryHCounters (const Molecule &query, Array<int> &counters, const int *mapping)
{
   int i;

   if (mapping == 0)
      counters.clear_resize(query.vertexEnd());

   if (mapping == 0)
      for (i = query.vertexBegin(); i < query.vertexEnd(); i = query.vertexNext(i))
         counters[i] = 0;
   else
      for (i = query.vertexBegin(); i < query.vertexEnd(); i = query.vertexNext(i))
         counters[mapping[i]] = 0;

   for (i = query.vertexBegin(); i < query.vertexEnd(); i = query.vertexNext(i))
   {
      const Vertex &vertex = query.getVertex(i);
      const Atom &atom = query.getAtom(i);

      if (atom.label == ELEM_H)
      {
         if (vertex.degree() > 1)
            throw Error("hydrogen atom has degree %d", vertex.degree());

         if (vertex.degree() > 0)
         {
            int nei_idx = vertex.neiVertex(vertex.neiBegin());

            counters[mapping == 0 ? nei_idx: mapping[nei_idx]]++;
         }
      }

      if (query.getStereocenters().getType(i) > 0)
      {
         if (vertex.degree() == 3)
         {
            if ((atom.label == ELEM_C  && atom.charge == 0 && atom.radical == 0) ||
                (atom.label == ELEM_Si && atom.charge == 0 && atom.radical == 0) ||
                (atom.label == ELEM_N  && atom.charge == 1 && atom.radical == 0))
               counters[mapping == 0 ? i : mapping[i]]++;
         }
      }
   }

   for (i = query.edgeBegin(); i < query.edgeEnd(); i = query.edgeNext(i))
   {
      const Edge &edge = query.getEdge(i);
      const QueryBond &query_bond = query.getQueryBond(i);

      if (query_bond.stereo_care)
      {
         const Atom &beg = query.getAtom(edge.beg);
         const Atom &end = query.getAtom(edge.end);

         bool beg_lost_h = ((beg.label == ELEM_C || beg.label == ELEM_Si) && beg.charge == 0 && beg.radical == 0) ||
                            (beg.label == ELEM_N && beg.charge == 1 && beg.radical == 0);
         bool end_lost_h = ((end.label == ELEM_C || end.label == ELEM_Si) && end.charge == 0 && end.radical == 0) ||
                            (end.label == ELEM_N && end.charge == 1 && end.radical == 0);

         if (query.getVertex(edge.beg).degree() < 3 && beg_lost_h)
            counters[mapping == 0 ? edge.beg : mapping[edge.beg]]++;
         if (query.getVertex(edge.end).degree() < 3 && end_lost_h)
            counters[mapping == 0 ? edge.end : mapping[edge.end]]++;
      }
   }
}

void MoleculeSubstructureMatcher::calcTargetHCounters (const Molecule &target, Array<int> &counters)
{
   int i;

   counters.clear_resize(target.vertexEnd());

   for (i = target.vertexBegin(); i < target.vertexEnd(); i = target.vertexNext(i))
   {
      counters[i] = target.getAtom(i).implicit_h;
      
      if (counters[i] < 0)
         counters[i] = 0;
   }

   for (i = target.vertexBegin(); i < target.vertexEnd(); i = target.vertexNext(i))
   {
      const Vertex &vertex = target.getVertex(i);
      const Atom &atom = target.getAtom(i);

      if (atom.label == ELEM_H)
      {
         if (vertex.degree() > 1)
            throw Error("hydrogen atom has degree %d", vertex.degree());

         if (vertex.degree() > 0)
         {
            int nei_idx = vertex.neiVertex(vertex.neiBegin());

            counters[nei_idx]++;
         }
      }
   }
}


void MoleculeSubstructureMatcher::setQuery (const Molecule &query)
{
   int i;

   if (query.getRGroups() != 0 && query.getRGroups()->getRGroupCount() > 0)
   {
      _markush.reset(new MarkushContext(query, _target));
      _query = &_markush->query;
   } else
   {
      _markush.reset(0);
      _query = &query;
   }

   calcQueryHCounters(*_query, _query_h_counters, 0);
   
   QS_DEF(Array<int>, ignored);

   ignored.clear_resize(_query->vertexEnd());

   markIgnoredHydrogens(*_query, ignored.ptr(), 0, 1);

   _3d_constrained_atoms.clear_resize(_query->vertexEnd());
   _3d_constrained_atoms.zerofill();

   if (_query->get3dConstraints() != 0)
   {
      Molecule3dConstraintsChecker checker(*_query->get3dConstraints());

      checker.markUsedAtoms(_3d_constrained_atoms.ptr(), 1);
   }

   if (_shouldUnfoldTargetHydrogens(*_query))
   {
      _did_h_unfold = true;
      _target.unfoldHydrogens(&_unfolded_target_h);
      calcTargetHCounters(_target, _target_h_counters);
   }
   else
      _did_h_unfold = false;

   _ee.create(_target);
   _ee->cb_match_vertex = _matchAtoms;
   _ee->cb_match_edge = _matchBonds;
   _ee->cb_vertex_add = _addAtom;
   _ee->cb_vertex_remove = _removeAtom;
   _ee->cb_edge_add = _addBond;
   _ee->cb_embedding = _embedding;
   _ee->userdata = this;

   _ee->setSubgraph(*_query);
   for (i = _query->vertexBegin(); i != _query->vertexEnd(); i = _query->vertexNext(i))
   {
      if ((ignored[i] && !_3d_constrained_atoms[i]) ||
          (_query->haveQueryAtoms() && _query->getQueryAtom(i).type == QUERY_ATOM_RGROUP))
         _ee->ignoreSubgraphVertex(i);
   }


}

void MoleculeSubstructureMatcher::setNeiCounters (
   const MoleculeAtomNeighbourhoodCounters *query_counters, 
   const MoleculeAtomNeighbourhoodCounters *target_counters)
{
   _query_nei_counters = query_counters;
   _target_nei_counters = target_counters;
}

bool MoleculeSubstructureMatcher::find ()
{
   if (_query == 0)
      throw Error("no query");

   if (match_3d != 0 && !_query->have_xyz)
      throw Error("cannot do 3D match without XYZ in the query");

   if (match_3d != 0 && !_target.have_xyz)
      return false;

   if (highlighting != 0)
      highlighting->init(_target);

   if (_canUseEquivalenceHeuristic(*_query))
      _ee->setEquivalenceHandler(vertex_equivalence_handler);
   else
      _ee->setEquivalenceHandler(NULL);

   Obj<AromaticityMatcher> am;
   if (use_aromaticity_matcher)
   {
      am.create(*_query, _target);
      _am = am.get();
   }
   else
      _am = 0;

   Obj<MoleculePiSystemsMatcher> pi_systems_matcher;
   if (use_pi_systems_matcher)
   {
      pi_systems_matcher.create(_target);
      _pi_systems_matcher = pi_systems_matcher.get();
   }
   else
      _pi_systems_matcher = 0;

   int result = _ee->process();

   if (_did_h_unfold)
   {
      QS_DEF(Array<int>, atoms_to_remove);
      atoms_to_remove.clear();
      for (int i = 0; i < _unfolded_target_h.size(); i++)
      {
         if (_unfolded_target_h[i])
         {
            // Update implicit hygrogens count
            const Vertex &v = _target.getVertex(i);
            int nei_index = v.neiVertex(v.neiBegin());

            _target.getAtom2(nei_index).implicit_h++;

            // Remove from highlighting
            if (highlighting != 0)
               highlighting->removeVertex(i);
            atoms_to_remove.push(i);
         }
      }
      if (atoms_to_remove.size() > 0)
         _target.removeAtoms(atoms_to_remove);
   }

   return result == 0;
}

int MoleculeSubstructureMatcher::_countSubstituents (const Molecule &mol, int idx)
{
   int i, res = 0;

   const Vertex &vertex = mol.getVertex(idx);

   for (i = vertex.neiBegin(); i != vertex.neiEnd(); i = vertex.neiNext(i))
   {
      const Atom &nei_atom = mol.getAtom(vertex.neiVertex(i));

      if (nei_atom.label != ELEM_H)
         res++;
   }

   return res;
}

bool MoleculeSubstructureMatcher::matchAtoms (const Graph &subgraph, const Graph &supergraph,
                                               const int *core_sub, int sub_idx, int super_idx, 
                                               dword flags)
{
   const Molecule &query = (const Molecule &)subgraph;
   const Molecule &target  = (const Molecule &)supergraph;

   const Vertex &sub_vertex = query.getVertex(sub_idx);
   const Vertex &super_vertex = target.getVertex(super_idx);
   const Atom &super_atom = target.getAtom(super_idx);
   const Atom &sub_atom = query.getAtom(sub_idx);
   const QueryAtom &sub_query_atom = query.getQueryAtom(sub_idx);

   switch (sub_query_atom.type)
   {
      case QUERY_ATOM_A:
         if (super_atom.label == ELEM_H)
            return false;
         break;
      case QUERY_ATOM_Q:
         if (super_atom.label == ELEM_H || super_atom.label == ELEM_C)
            return false;
         break;
      case QUERY_ATOM_X:
         if (!Element::isHalogen(super_atom.label))
            return false;
         break;
      case QUERY_ATOM_LIST:
         if (!sub_query_atom.listContains(super_atom.label))
            return false;
         break;
      case QUERY_ATOM_NOTLIST:
         if (sub_query_atom.listContains(super_atom.label))
            return false;
         break;
      default:
      {
         if (query.isPseudoAtom(sub_idx) && target.isPseudoAtom(super_idx))
         {
            if (strcmp(query.getPseudoAtom(sub_idx), target.getPseudoAtom(super_idx)) != 0)
               return false;
         }
         else if (!query.isPseudoAtom(sub_idx) && !target.isPseudoAtom(super_idx))
         {
            if (sub_atom.label != super_atom.label)
               return false;
         }
         else
            return false;
      }
   }

   if (sub_atom.isotope != 0 && sub_atom.isotope != super_atom.isotope)
      return false;
   if (flags & MATCH_ATOM_CHARGE)
      if (sub_atom.charge != 0 && sub_atom.charge != super_atom.charge)
         return false;
   if (sub_atom.radical != 0 && sub_atom.radical != super_atom.radical)
      return false;
   if (flags & MATCH_ATOM_VALENCE)
      if (sub_atom.explicit_valence && sub_atom.valence != super_atom.valence)
         return false;

   if (query.getStereocenters().getType(sub_idx) > target.getStereocenters().getType(super_idx))
      return false;

   if (sub_query_atom.substitution_count != 0)
   {
      int subst = _countSubstituents(target, super_idx);

      if (sub_query_atom.substitution_count == -1)
      {
         if (subst != 0)
            return false;
      }
      else if (sub_query_atom.substitution_count == -2)
      {
         if (subst != _countSubstituents(query, sub_idx))
            return false;
      }
      else if (sub_query_atom.substitution_count < 6)
      {
         if (subst != sub_query_atom.substitution_count)
            return false;
      }
      else
      {
         if (subst < 6)
            return false;
      }
   }

   if (sub_query_atom.unsaturation)
   {
      int i;
      
      for (i = super_vertex.neiBegin(); i != super_vertex.neiEnd(); i = super_vertex.neiNext(i))
      {
         int edge_idx = super_vertex.neiEdge(i);

         if (target.getBond(edge_idx).type != BOND_SINGLE)
            break;
      }
      
      if (i == super_vertex.neiEnd())
         return false;
   }

   if (sub_query_atom.ring_bond_count != 0)
   {
      int i, rings = 0;

      for (i = super_vertex.neiBegin(); i != super_vertex.neiEnd(); i = super_vertex.neiNext(i))
         if (target.bondInRing(super_vertex.neiEdge(i)))
            rings++;
      
      if (sub_query_atom.ring_bond_count == -1)
      {
         if (rings != 0)
            return false;
      }
      else if (sub_query_atom.ring_bond_count == -2)
      {
         int rings_sub = 0;

         for (i = sub_vertex.neiBegin(); i != sub_vertex.neiEnd(); i = sub_vertex.neiNext(i))
         {
            if (query.bondInRing(sub_vertex.neiEdge(i)))
               rings_sub++;
         }
         if (rings != rings_sub)
            return false;
      }
      else if (sub_query_atom.ring_bond_count < 4)
      {
         if (rings != sub_query_atom.ring_bond_count)
            return false;
      }
      else
      {
         if (rings < 4)
            return false;
      }
   }

   if (sub_query_atom.aromaticity != 0)
   {
      bool arom = false;
      int i;

      for (i = super_vertex.neiBegin(); i != super_vertex.neiEnd(); i = super_vertex.neiNext(i))
         if (target.getBond(super_vertex.neiEdge(i)).type == BOND_AROMATIC)
         {
            arom = true;
            break;
         }

      if ((sub_query_atom.aromaticity == QUERY_ATOM_AROMATIC) && !arom)
         return false;

      if ((sub_query_atom.aromaticity == QUERY_ATOM_ALIPHATIC) && arom)
         return false;

   }

   return true;
}

bool MoleculeSubstructureMatcher::_matchAtoms(const Graph &subgraph, const Graph &supergraph,
                                              const int *core_sub, int sub_idx, int super_idx, void *userdata)
{
   MoleculeSubstructureMatcher *self = (MoleculeSubstructureMatcher *)userdata;

   if (self->_did_h_unfold && (&subgraph == (const Graph *)self->_query))
   {
      if (sub_idx < self->_3d_constrained_atoms.size() && self->_3d_constrained_atoms[sub_idx])
         // we can't check 3D constraint on unfolded atom, because it has no actual position
         if (self->_unfolded_target_h[super_idx])
            return false;
   }

   dword match_atoms_flags = 0xFFFFFFFF;
   // If target atom belongs to a pi-system then its charge
   // should be checked after embedding
   if (self->_pi_systems_matcher)
   {
      if (self->_pi_systems_matcher->isAtomInPiSystem(super_idx))
         match_atoms_flags &= ~(MATCH_ATOM_CHARGE | MATCH_ATOM_VALENCE);
   }

   if (!matchAtoms(subgraph, supergraph, core_sub, sub_idx, super_idx, match_atoms_flags))
      return false;

   if (sub_idx < self->_query_h_counters.size() && super_idx < self->_target_h_counters.size())
      if (self->_query_h_counters[sub_idx] > self->_target_h_counters[super_idx])
         return false;

   const QueryAtom &sub_query_atom = ((const Molecule &)subgraph).getQueryAtom(sub_idx);

   if (sub_query_atom.hydrogens == 1) // means H0
   {
      if (self->_target_h_counters[super_idx] > self->_query_h_counters[sub_idx])
         return false;
   }
   else if (sub_query_atom.hydrogens > 1) // means Hn
   {
      if (self->_target_h_counters[super_idx] < (sub_query_atom.hydrogens - 1))
         return false;
   }

   if (self->_query_nei_counters != 0 && self->_target_nei_counters != 0)
   {
      bool use_bond_types = (self->_pi_systems_matcher == 0);
      bool ret = self->_query_nei_counters->testSubstructure(
         *self->_target_nei_counters, sub_idx, super_idx, use_bond_types);
      if (!ret)
         return false;
   }

   if (self->match_3d == AFFINE)
   {
      QS_DEF(Array<int>, core_sub_full);

      core_sub_full.copy(core_sub, subgraph.vertexEnd());
      core_sub_full[sub_idx] = super_idx;

      GraphAffineMatcher matcher(subgraph, supergraph, core_sub_full.ptr());

      matcher.cb_get_xyz = getAtomPos;

      const Molecule &query = (const Molecule &)subgraph;
      int total_fixed = query.vertexCount();

      if (query.fixed_atoms.size() > 0)
      {
         matcher.fixed_vertices = &query.fixed_atoms;
         total_fixed = query.fixed_atoms.size();
      }

      if (!matcher.match(self->rms_threshold * sqrt((float)total_fixed)))
         return false;
   }

   return true;
}

bool MoleculeSubstructureMatcher::matchBonds (const Graph &subgraph, const Graph &supergraph,
                                               int sub_idx, int super_idx, AromaticityMatcher *am,
                                               dword flags)
{
   const Molecule &query = (const Molecule &)subgraph;
   const Molecule &molecule = (const Molecule &)supergraph;

   const Bond &super_bond = molecule.getBond(super_idx);
   const Bond &sub_bond = query.getBond(sub_idx);
   const QueryBond &sub_query_bond = query.getQueryBond(sub_idx);

   if (flags & MATCH_BOND_TYPE)
   {
      if (sub_query_bond.type == 0) // ordinary (not query) bond type;
      {
         if (super_bond.type != BOND_AROMATIC)
         {
            if (super_bond.type != sub_bond.type)
               return false;
         }
         else if (sub_bond.type != BOND_AROMATIC)
         {
            if (am == 0)
               return false;

            // It's incorrect to test query bond in ring, because of markush search:
            //    cycle part can be in another R-Group fragment
            if (!am->isAbleToFixBond(sub_idx))
               return false;
         }
      }
      else switch (sub_query_bond.type)
      {
         case QUERY_BOND_SINGLE_OR_DOUBLE:
            if (super_bond.type != BOND_SINGLE && super_bond.type != BOND_DOUBLE)
               return false;
            break;
         case QUERY_BOND_SINGLE_OR_AROMATIC:
            if (super_bond.type != BOND_SINGLE && super_bond.type != BOND_AROMATIC)
               return false;
            break;
         case QUERY_BOND_DOUBLE_OR_AROMATIC:
            if (super_bond.type != BOND_DOUBLE && super_bond.type != BOND_AROMATIC)
               return false;
            break;
      }
   }

   if (sub_query_bond.topology != 0)
   {
      bool in_ring = molecule.bondInRing(super_idx);

      if (sub_query_bond.topology == QUERY_BOND_TOPOLOGY_CHAIN && in_ring)
         return false;
      if (sub_query_bond.topology == QUERY_BOND_TOPOLOGY_RING && !in_ring)
         return false;
   }

   return true;
}

bool MoleculeSubstructureMatcher::_matchBonds (const Graph &subgraph, const Graph &supergraph,
                                               int sub_idx, int super_idx, void *userdata)
{
   MoleculeSubstructureMatcher *self = (MoleculeSubstructureMatcher *)userdata;

   dword flags = 0xFFFFFFFF;
   // If target bond belongs to a pi-system then it
   // should be checked after embedding
   if (self->_pi_systems_matcher)
   {
      if (self->_pi_systems_matcher->isBondInPiSystem(super_idx))
         flags &= ~MATCH_BOND_TYPE;
   }

   return matchBonds(subgraph, supergraph, sub_idx, super_idx, self->_am, flags);
}

void MoleculeSubstructureMatcher::removeAtom (const Graph &subgraph, int sub_idx, AromaticityMatcher *am)
{
   if (am == 0)
      return;

   am->unfixNeighbourAromaticBonds(sub_idx);
}

void MoleculeSubstructureMatcher::_removeAtom (const Graph &subgraph, int sub_idx, void *userdata)
{
   MoleculeSubstructureMatcher *self = (MoleculeSubstructureMatcher *)userdata;
   
   removeAtom(subgraph, sub_idx, self->_am);
}

void MoleculeSubstructureMatcher::addAtom (const Graph &subgraph, const Graph &supergraph,
                                            int sub_idx, int super_idx, AromaticityMatcher *am)
{
   if (am == 0)
      return;

   am->fixAtom(sub_idx, super_idx);
}

void MoleculeSubstructureMatcher::_addAtom (const Graph &subgraph, const Graph &supergraph,
                                            int sub_idx, int super_idx, void *userdata)
{
   MoleculeSubstructureMatcher *self = (MoleculeSubstructureMatcher *)userdata;

   addAtom(subgraph, supergraph, sub_idx, super_idx, self->_am);
}

void MoleculeSubstructureMatcher::addBond (const Graph &subgraph, const Graph &supergraph,
                                            int sub_idx, int super_idx, AromaticityMatcher *am)
{
   if (am == 0)
      return;

   // Check bonds types
   const Molecule &query = (const Molecule &)subgraph;
   const Molecule &molecule = (const Molecule &)supergraph;

   const Bond &super_bond = molecule.getBond(super_idx);
   const Bond &sub_bond = query.getBond(sub_idx);

   if (super_bond.type == BOND_AROMATIC && sub_bond.type != BOND_AROMATIC)
      am->fixAromaticBond(sub_idx);
}

void MoleculeSubstructureMatcher::_addBond (const Graph &subgraph, const Graph &supergraph,
                                            int sub_idx, int super_idx, void *userdata)
{
   MoleculeSubstructureMatcher *self = (MoleculeSubstructureMatcher *)userdata;

   addBond(subgraph, supergraph, sub_idx, super_idx, self->_am);
}

int MoleculeSubstructureMatcher::_embedding (const Graph &subgraph, const Graph &supergraph,
                                             int *core_sub, int *core_super, void *userdata)
{
   MoleculeSubstructureMatcher *self = (MoleculeSubstructureMatcher *)userdata;

   if (self->_markush.get() == 0 || self->_markush->sites.size() == self->_markush->depth)
      return self->_embedding_common(core_sub, core_super);
   else
      return self->_embedding_markush(core_sub, core_super);
}

int MoleculeSubstructureMatcher::_embedding_common (int *core_sub, int *core_super)
{
   const Molecule &query = *_query;

   if (!MoleculeStereocenters::checkSub(query.getStereocenters(),
                                       _target.getStereocenters(), core_sub, false))
      return 1;

   if (!MoleculeCisTrans::checkSub(query, _target, core_sub))
      return 1;

   if (query.get3dConstraints() != 0)
   {
      Molecule3dConstraintsChecker checker(*query.get3dConstraints());

      if (!checker.check(_target, core_sub))
         return 1;
   }

   // Check possible aromatic configuration
   if (_am != 0)
   {
      _am->setQueryVertexIgnoreFilterPtr(core_sub, EmbeddingEnumerator::IGNORE);
      if (!_am->isAbleToMatch())
         return 1;
   }

   // Check possible pi-systems configurations
   if (_pi_systems_matcher != 0)
   {
      if (!_pi_systems_matcher->checkEmbedding(query, core_sub))
         return 1;
   }

   // affine transformation match
   /*if (match_3d == AFFINE)
   {
      GraphAffineMatcher matcher(query, _target, core_sub);

      matcher.cb_get_xyz = getAtomPos;

      if (query.fixed_atoms.size() > 0)
         matcher.fixed_vertices = &query.fixed_atoms;

      if (!matcher.match(v1, v2))
         return 1;
   }*/

   // chemical conformation match
   if (match_3d == CONFORMATION)
   {
      QS_DEF(Array<int>, mapping);
      GraphDecomposer decomposer(query);

      Filter filter(core_sub, Filter::MORE, -1);
      decomposer.decompose(&filter);

      int i, comp;

      for (comp = 0; comp < decomposer.getComponentsCount(); comp++)
      {
         mapping.clear_resize(query.vertexEnd());
         mapping.fffill();
         for (i = query.vertexBegin(); i != query.vertexEnd(); i = query.vertexNext(i))
         {
            if (decomposer.getComponent(i) == comp)
               mapping[i] = core_sub[i];
         }

         EdgeRotationMatcher matcher(query, _target, mapping.ptr());

         matcher.cb_get_xyz = getAtomPos;
         matcher.cb_can_rotate = _isSingleBond;
         matcher.equalize_edges = true;

         if (!matcher.match(rms_threshold, 0.1f))
            return 1;
      }
   }

   if (_markush.get() != 0)
      if (!_checkRGroupConditions())
         return 1;

   if (highlighting != 0)
      highlighting->onSubgraph(query, core_sub);
   
   if (find_all_embeddings)
   {
      if (cb_embedding != 0)
         cb_embedding(query, _target, core_sub, core_super, cb_embedding_context);

      return 1;
   }

   return 0;
}

int MoleculeSubstructureMatcher::_embedding_markush (int *core_sub, int *core_super)
{
   const Molecule &g1 = *_query;
   MarkushContext &context = *_markush.get();

   // Current site atom index
   int old_site_idx = context.sites[context.depth];
   bool two_att_points;

   // Check number of attachment points for current site
   switch (g1.getVertex(old_site_idx).degree())
   {
   case 1:
      two_att_points = false;
      break;
   case 2:
      two_att_points = true;
      break;
   default:
      return true;
   }

   // For all possible rgroups at current site
   // do not do this:
   // const Array<int> &rg_list = g1.getRGroups()->getSiteRGroups(old_site_idx);

   for (int rg_idx = 0; rg_idx < g1.getRGroups()->getSiteRGroups(old_site_idx).size(); rg_idx++)
   {
      const RGroup &rgroup = g1.getRGroups()->getRGroup(g1.getRGroups()->getSiteRGroups(old_site_idx)[rg_idx]);
      // For all rgroup fragments
      for (int fr_idx = 0; fr_idx < rgroup.fragments.size(); fr_idx++)
      {
         const Molecule &fragment = *rgroup.fragments[fr_idx];

         if (fragment.getRGroupFragment()->attachmentPointCount() > 2)
            throw Error("more than two attachment points");

         MoleculeRGroupFragment &att_fragment = *fragment.getRGroupFragment();

         int att_idx1;
         int att_idx2;
         int i, j;

         for (i = 0; (att_idx1 = att_fragment.getAttachmentPoint(0, i)) != -1; i++)
         {
            if (two_att_points)
            {
               for (j = 0; (att_idx2 = att_fragment.getAttachmentPoint(1, j)) != -1; j++)
                  if (!_attachRGroupAndContinue(core_sub, core_super, &fragment, true, att_idx1, att_idx2, g1.getRGroups()->getSiteRGroups(old_site_idx)[rg_idx]))
                     return 0;
            } else if (!_attachRGroupAndContinue(core_sub, core_super, &fragment, false, att_idx1, -1, g1.getRGroups()->getSiteRGroups(old_site_idx)[rg_idx]))
               return 0;
         }
      }
   }

   if (!two_att_points && !_attachRGroupAndContinue(core_sub, core_super, 0, false, -1, -1, -1))
      return 0;

   return 1;
}

bool MoleculeSubstructureMatcher::_attachRGroupAndContinue (int *core1, int *core2,
                                               const Molecule *fragment, bool two_attachment_points,
                                               int att_idx1, int att_idx2, int rgroup_idx)
{
   MarkushContext &context = *_markush.get();
   QS_DEF(Array<int>, fr_mapping);
   int i;
   int n_sites = 0;
   bool ok = true;
   int hydrogen_attached_to = -1;

   // Site atom index in new molecule
   int cur_site = context.sites[context.depth];

   int src_att_idx1 = -1, src_att_idx2 = -1;
   Bond rg_bond1(0), rg_bond2(0);
   QueryBond rg_qbond1(0), rg_qbond2(0);

   // Paramaters for stereocenter restoration
   bool stereo_was_saved = false;
   int saved_stereo_type, saved_stereo_group, saved_stereo_pyramid[4];

   // If some rgroup must be attached
   if (rgroup_idx != -1)
   {
      // Add rgroup fragment to new molecule
      context.query.mergeWithMolecule(*fragment, &fr_mapping);

      _query_h_counters.expand(context.query.vertexEnd());

      calcQueryHCounters(*fragment, _query_h_counters, fr_mapping.ptr());

      // Site atom vertex in new molecule
      const Vertex &cur_site_vertex = context.query.getVertex(cur_site);

      // First and second neighbors indices and neighboring atoms indices
      int nei_idx1 = cur_site_vertex.neiBegin();
      int nei_idx2 = cur_site_vertex.neiNext(nei_idx1);

      int nei_atom_idx1 = context.query.getRGroups()->getAttachmentOrder(cur_site, 0);
      int nei_atom_idx2 = context.query.getRGroups()->getAttachmentOrder(cur_site, 1);

      if (two_attachment_points)
      {
         if (cur_site_vertex.neiVertex(nei_idx1) != nei_atom_idx1)
         {
            // Swap attachment neighbors according to attachment orders (first and second)
            __swap(nei_idx1, nei_idx2, i);
         }

         if (cur_site_vertex.neiVertex(nei_idx1) != nei_atom_idx1 ||
            cur_site_vertex.neiVertex(nei_idx2) != nei_atom_idx2)
            throw Error("Internal or rgroup query file error!");
      } else {
         if (att_idx1 == -1)
         {
            // Second attachment index of fragment must be attached to second-ordered neighbor
            if (cur_site_vertex.neiVertex(nei_idx1) != nei_atom_idx2)
               ok = false;
            att_idx1 = att_idx2;
         } else if (cur_site_vertex.neiVertex(nei_idx1) != nei_atom_idx1)
            // First attachment index of fragment must be attached to first-ordered neighbor
            ok = false;
      }

      if (ok)
      {
         // First attachment bond of current rgroup site
         src_att_idx1 = cur_site_vertex.neiVertex(nei_idx1);
         context.query.flipBond(src_att_idx1, cur_site, fr_mapping[att_idx1]);

         if (two_attachment_points)
         {
            // Second attachment bond of current rgroup site
            src_att_idx2 = cur_site_vertex.neiVertex(nei_idx2);
            context.query.flipBond(src_att_idx2, cur_site, fr_mapping[att_idx2]);
         } 
         else
         {
            // Check if there is only explicit hydrogen attaching
            if (fragment->vertexCount() == 1 &&
               fragment->getAtom(fragment->vertexBegin()).label == ELEM_H)
            {
               // Increase explicit hydrogens count
               hydrogen_attached_to = src_att_idx1;
               _query_h_counters[hydrogen_attached_to]++;

               // Check if the atom with changed hydrogens count is still matching
               if (!_matchAtoms(context.query, _target, core1, src_att_idx1, core1[src_att_idx1], this))
                  ok = false;
            }
         }
      }
   }
   else 
   {                   
      // no rgroup (rgroup = implicit hydrogen)
      // This might happen only for rgroup with one attachment point
      const Vertex &cur_site_vertex = context.query.getVertex(cur_site);

      if (cur_site_vertex.degree() > 0 && _target.getAtom(core1[cur_site_vertex.neiVertex(cur_site_vertex.neiBegin())]).implicit_h == 0)
         return true;

      int nei_idx1 = cur_site_vertex.neiBegin();
      src_att_idx1 = cur_site_vertex.neiVertex(nei_idx1);

      // Remove edge to site
      rg_qbond1 = context.query.getQueryBond(cur_site_vertex.neiEdge(nei_idx1));
      rg_bond1 = context.query.getBond(cur_site_vertex.neiEdge(nei_idx1));

      // Save stereocenter before bond removing because bond 
      // removing can cause stereocenter destruction
      MoleculeStereocenters &qstereo = context.query.getStereocenters();
      if (qstereo.exists(src_att_idx1))
      {
         qstereo.get(src_att_idx1, saved_stereo_type, saved_stereo_group, saved_stereo_pyramid);
         stereo_was_saved = true;
      }

      context.query.removeBond(cur_site_vertex.neiEdge(nei_idx1));
   }

   // Init new embedding enumerator context
   EmbeddingEnumerator ee(_target);

   ee.cb_match_vertex = _matchAtoms;
   ee.cb_match_edge = _matchBonds;
   ee.cb_vertex_add = _addAtom;
   ee.cb_vertex_remove = _removeAtom;
   ee.cb_edge_add = _addBond;
   ee.cb_embedding = _embedding;

   ee.setSubgraph(context.query);

   ee.userdata = this;

   if (rgroup_idx != -1)
   {
      context.query_marking.expand(context.query.vertexEnd());

      // Check if there are explicit hydrogens or new sites attached- ignore them in further matching
      // Mark new vertices and add new sites if any
      QS_DEF(Array<int>, ignored);

      ignored.clear_resize(fragment->vertexEnd());

      markIgnoredHydrogens(*fragment, ignored.ptr(), 0, 1);

      if (hydrogen_attached_to != -1)
      {
         int idx = fragment->vertexBegin();

         context.query_marking[fr_mapping[idx]] = context.depth;

         if (fragment->getAtom(idx).pureHydrogen())
            ee.ignoreSubgraphVertex(fr_mapping[idx]);
      } else
      {
         for (i = fragment->vertexBegin(); i < fragment->vertexEnd(); i = fragment->vertexNext(i))
         {
            context.query_marking[fr_mapping[i]] = context.depth;

            if (ignored[i])
               ee.ignoreSubgraphVertex(fr_mapping[i]);
            else if (fragment->getQueryAtom(i).type == QUERY_ATOM_RGROUP)
            {
               context.sites.push(fr_mapping[i]);
               n_sites++;
               ee.ignoreSubgraphVertex(fr_mapping[i]);
            }
         }
      }

      if (_am != 0)
         _am->validateQuery();
   }

   // Copy ingnored query vertices from previous state
   for (i = context.query.vertexBegin(); i < context.query.vertexEnd(); i = context.query.vertexNext(i))
   {
      if (context.query_marking[i] == context.depth)
         continue;

      if (core1[i] == EmbeddingEnumerator::IGNORE)
         ee.ignoreSubgraphVertex(i);
   }

   // Copy ignored target vertices
   for (i = _target.vertexBegin(); i < _target.vertexEnd(); i = _target.vertexNext(i))
      if (core2[i] == EmbeddingEnumerator::IGNORE)
         ee.ignoreSupergraphVertex(i);

   // Copy partial mapping from previous state
   for (i = context.query.vertexBegin(); i < context.query.vertexEnd(); i = context.query.vertexNext(i))
   {
      if (context.query_marking[i] == context.depth)
         continue;

      if (core1[i] >= 0 && !ee.fix(i, core1[i]))
         ok = false;
   }

   // Store attached rgroup index at current site
   context.sites[context.depth] = rgroup_idx;

   context.depth++;

   // Continue in new state
   if (ok && !ee.process())
      return false;

   // Restore state
   context.depth--;

   context.sites[context.depth] = cur_site;

   if (rgroup_idx != -1)
   {
      context.query.flipBond(src_att_idx1, fr_mapping[att_idx1], cur_site);

      if (src_att_idx2 != -1)
         context.query.flipBond(src_att_idx2, fr_mapping[att_idx2], cur_site);
   } 
   else // no rgroup
   {
      // Restore edge to site
      int bond_added = context.query.addBondCopy(src_att_idx1, cur_site, rg_bond1);
      if (context.query.haveQueryBonds())
         context.query.getQueryBond2(bond_added) = rg_qbond1;

      // Restore stereocenter
      if (stereo_was_saved)
      {
         MoleculeStereocenters &qstereo = context.query.getStereocenters();
         if (qstereo.exists(src_att_idx1))
            qstereo.remove(src_att_idx1);
         qstereo.add(src_att_idx1, saved_stereo_type, saved_stereo_group, saved_stereo_pyramid);
      }
   }

   if (rgroup_idx != -1)
   {
      Filter remove_filter(context.query_marking.ptr(), Filter::EQ, context.depth);
      context.query.removeAtoms(remove_filter);

      if (hydrogen_attached_to != -1)
         _query_h_counters[hydrogen_attached_to]--;

      while (n_sites-- > 0)
         context.sites.pop();

      if (_am != 0)
         _am->validateQuery();
   }

   return true;
}

bool MoleculeSubstructureMatcher::_checkRGroupConditions ()
{
   QS_DEF(Array<int>, occurences);
   QS_DEF(Array<int>, conditions);
   QS_DEF(Array<int>, nontop_rgroups);
   bool resth_satisfied;
   int i, k;
   MarkushContext &context = *_markush.get();
   MoleculeRGroups &rgroups = *_query->getRGroups();

   int n_rgroups = rgroups.getRGroupCount();

   occurences.clear_resize(n_rgroups);
   occurences.zerofill();
   conditions.clear_resize(n_rgroups);
   conditions.zerofill();
   nontop_rgroups.clear_resize(n_rgroups);
   nontop_rgroups.zerofill();

   // Count occurences of each rgroup
   for (i = 0; i < context.sites.size(); i++)
      if (context.sites[i] >= 0)
         occurences[context.sites[i]]++;


   for (int i = 0; i < n_rgroups; i++)
   {
      const RGroup & rgroup = rgroups.getRGroup(i);

      if (rgroup.fragments.size() == 0)
      {
         conditions[i] = -1;
         continue;
      }

      if (rgroup.if_then >= 0)
         nontop_rgroups[rgroup.if_then] = 1;

      // Check occurence
      if (!rgroup.occurenceSatisfied(occurences[i]))
         continue;

      resth_satisfied = true;

      // Check RestH
      if (rgroup.rest_h == 1)
      {
         for (int site_idx = rgroups.begin(), j = 0; site_idx < rgroups.end() && resth_satisfied; site_idx = rgroups.next(site_idx), j++)
         {
            const Array<int> &rg_list = rgroups.getSiteRGroups(rgroups.atomIdx(site_idx));

            for (k = 0; k < rg_list.size() && resth_satisfied; k++)
               if (rg_list[k] == i)
               {
                  if (context.sites[j] != i)
                     resth_satisfied = false;
                  break;
               }
         }
         if (!resth_satisfied)
            continue;
      }

      // Without if/then
      conditions[i] = 1;
   }

   // Check if/then relations
   for (int i = 0; i < n_rgroups; i++)
   {
      if (nontop_rgroups[i] == 1 || conditions[i] == -1)
         continue;

      const RGroup &rgroup = rgroups.getRGroup(i);
      const RGroup *iter = &rgroup;
      int i_iter = i;
      int counter = 0;

      do 
      {
         if (conditions[i_iter] != 1 || counter >= n_rgroups)
            return false;

         if ((i_iter = iter->if_then) >= 0)
            iter = &rgroups.getRGroup(i_iter);

         counter++;
      } while(i_iter >= 0);
   }

   return true;
}

const int * MoleculeSubstructureMatcher::getQueryMapping ()
{
   return _ee->getSubgraphMapping();
}

const int * MoleculeSubstructureMatcher::getTargetMapping ()
{
   return _ee->getSupergraphMapping();
}

void MoleculeSubstructureMatcher::ignoreQueryAtom (int idx)
{
   _ee->ignoreSubgraphVertex(idx);
}

void MoleculeSubstructureMatcher::ignoreTargetAtom (int idx)
{
   _ee->ignoreSupergraphVertex(idx);
}

bool MoleculeSubstructureMatcher::fix (int query_atom_idx, int target_atom_idx)
{
   return _ee->fix(query_atom_idx, target_atom_idx);
}

void MoleculeSubstructureMatcher::markIgnoredHydrogens (const Molecule &mol, int *arr, int value_keep, int value_ignore)
{
   int i;

   for (i = mol.vertexBegin(); i != mol.vertexEnd(); i = mol.vertexNext(i))
      arr[i] = value_keep;

   for (i = mol.vertexBegin(); i != mol.vertexEnd(); i = mol.vertexNext(i))
   {
      if (!mol.getAtom(i).pureHydrogen())
         continue;

      const Vertex &vertex = mol.getVertex(i);

      if (vertex.degree() == 1)
      {
         int nei_idx = vertex.neiVertex(vertex.neiBegin());

         if (mol.getAtom(nei_idx).label != ELEM_H) // do not ignore rare H-H fragment
            arr[i] = value_ignore;
      }
   }

   // keep hydrogens that are parts of 3d constraints
   if (mol.get3dConstraints() != 0)
   {
      Molecule3dConstraintsChecker checker(*mol.get3dConstraints());

      checker.markUsedAtoms(arr, value_keep);
   }
}

void MoleculeSubstructureMatcher::getAtomPos (const Graph &graph, int vertex_idx, Vec3f &pos)
{
   pos.copy(((const Molecule &)graph).getAtom(vertex_idx).pos);
}

bool MoleculeSubstructureMatcher::_isSingleBond (const Graph &graph, int edge_idx)
{
   return ((const Molecule &)graph).getBond(edge_idx).type == BOND_SINGLE;
}

bool MoleculeSubstructureMatcher::needCoords (int match_3d, const Molecule &query)
{
   if (match_3d != 0)
      return true;

   if (query.get3dConstraints() != 0)
      return true;

   return false;
}

bool MoleculeSubstructureMatcher::_canUseEquivalenceHeuristic (const Molecule &query)
{
   if (query.get3dConstraints() != 0)
      return false;
   if (query.getStereocenters().size() > 0)
      return false;
   return true;
}

